# "Partisanship, political awareness, and retrospective evaluations, 1956-2016"
# Phil Jones
# 02.26.19

# Replication file to run multi-level regression models shown in Table 3, and simulated first differences shown in Figure 3. 

library(car)
library(lme4)
library(Zelig)
library(ZeligChoice)
library(ggplot2)
library(texreg)
library(MCMCglmm)
library(dplyr)
library(xtable)

load("pejones_PB_multi.RData")



########################################################################
#  Set-up before models
########################################################################

# load functions to derive first differences from MCMCglmm models
source("prediction_MCMCglmm.R")


# create clean versions of data to use in models
na.data <- subset(m.data, select=c("retro_eval", "age", "white", "income", "educ", "female", "pidi", "year", "retro_eval_group", "id2", "time.m", "economic", "pocketbook","polar_congress.m", "change_gdp.m", "change_gdp.m2", "pres_election", "unified_govt"))
na.data <- na.omit(na.data)

na.data.aware <- subset(m.data, select=c("retro_eval", "age", "white", "income", "educ", "female", "pidi", "year", "retro_eval_group", "id2", "time.m", "economic", "pocketbook","polar_congress.m", "change_gdp.m", "change_gdp.m2", "pres_election", "unified_govt", "awareness"))
na.data.aware <- na.omit(na.data.aware)

na.data.econ <- subset(na.data, retro_eval_group=="1Sociotropic"|retro_eval_group=="Pocketbook")



########################################################################
#  Models
########################################################################


prior1 <- list(R=list(V=1,fix=1), G=list(G1=list(V=diag(3), nu=.002), G2=list(V=1, nu=0.002)))

MC1 <- MCMCglmm(as.factor(retro_eval)~
# indiv-level vars
age+white+income+educ+female+pidi+
# eval level var
retro_eval_group+
# year-level vars
time.m+polar_congress.m+pres_election+change_gdp.m+change_gdp.m2+unified_govt+
# year-indiv interactions
time.m*pidi+polar_congress.m*pidi+pres_election*pidi+change_gdp.m*pidi+change_gdp.m2*pidi+unified_govt*pidi,
# grouping factors
random=~us(pidi):year+id2,
data=na.data,
prior=prior1,
nitt=50000, burnin=5000, thin=10,
family="ordinal", pr=T)
summary(MC1)



MC2 <- MCMCglmm(as.factor(retro_eval)~
# indiv-level vars
age+white+income+educ+female+pidi+
# eval level var
retro_eval_group+
# year-level vars
time.m+polar_congress.m+pres_election+change_gdp.m+change_gdp.m2+unified_govt+
# year-indiv interactions
time.m*pidi+polar_congress.m*pidi+pres_election*pidi+change_gdp.m*pidi+change_gdp.m2*pidi+unified_govt*pidi,
# grouping factors
random=~us(pidi):year+id2,
data= na.data.econ,
prior=prior1,
nitt=50000, burnin=5000, thin=10,
family="ordinal", pr=T)
summary(MC2)


MC3 <- MCMCglmm(as.factor(retro_eval)~
# indiv-level vars
age+white+income+educ+female+pidi+pidi*awareness+
# eval level var
retro_eval_group+
# year-level vars
time.m+polar_congress.m+pres_election+change_gdp.m+change_gdp.m2+unified_govt+
# year-indiv interactions
time.m*pidi+polar_congress.m*pidi+pres_election*pidi+change_gdp.m*pidi+change_gdp.m2*pidi+unified_govt*pidi+
time.m*pidi*awareness+polar_congress.m*pidi*awareness+pres_election*pidi*awareness+change_gdp.m*pidi*awareness+change_gdp.m2*pidi*awareness+unified_govt*pidi*awareness,
# grouping factors
random=~us(pidi):year+id2,
data=na.data.aware,
prior=prior1,
nitt=50000, burnin=5000, thin=10,
family="ordinal", pr=T)
summary(MC3)


########################################################################
#  Simulations
########################################################################

# Range of values to simulate for polarization and GDP (mean-centered and original)
p.range <- seq(from=(.5-0.6440357), to=(.75-0.6440357), by=.01)
p.range.orig <- seq(from=.5, to=.75, by=.01)
e.range <- seq(from=min(m.data$change_gdp.m, na.rm=T), to=max(m.data$change_gdp.m, na.rm=T), by=.3)
e2.range <- e.range^2
e.range.orig <- seq(from=min(m.data$change_gdp, na.rm=T), to=max(m.data$change_gdp, na.rm=T), by=.3)


# Distribution of survey years for rug plots on figures
p.dat <- as.data.frame(cbind(unique(m.data$polar_congress),1:27))
e.dat <- as.data.frame(cbind(unique(m.data$change_gdp),1:27))





# FIGURE 3A

myX <- matrix(0, ncol=nrow(fixef(MC1)), nrow=1)
colnames(myX) <- row.names(fixef(MC1))
myX <- as.data.frame(myX)
myX$"(Intercept)" <- 1
myX$age <- 4.52
myX$white <- 1
myX$income <- 3
myX$educ <- 3
myX$female <- 1
myX$pres_election <- 1

myX1 <- myX
myX1$pidiOutParty <- 1
myX1$"pidiOutParty:pres_election" <- 1

fd.econ.MC1 <- matrix(NA, nrow=length(e.range), ncol=3)

for(i in 1:length(e.range)){
	OutX <- myX1
	OutX$change_gdp.m <- e.range[i]
	OutX$change_gdp.m2 <- e2.range[i]
	OutX$"pidiOutParty:change_gdp.m" <- e.range[i]
	OutX$"pidiOutParty:change_gdp.m2" <- e2.range[i]
	OutX <- as.matrix(OutX)
	
	InX <- myX
	InX$change_gdp.m <- e.range[i]
	InX$change_gdp.m2 <- e2.range[i]
	InX <- as.matrix(InX)

	ans <- fd.MCMCglmm.3cat(MC1, InX, OutX, Z=NULL, use="all", type="response")
	fd.econ.MC1[i,] <- as.matrix(ests.MCMCglmm.3cat(ans)[3,])
	
	print(i/length(e.range))

}

fd.econ.MC1 <- cbind(e.range, e.range.orig, fd.econ.MC1)
fd.econ.MC1 <- as.data.frame(fd.econ.MC1)
names(fd.econ.MC1) <- c("changeGDP.m", "changeGDP", "fd", "lo", "up")

g1 <- ggplot(data= fd.econ.MC1, aes(y=fd, x= changeGDP, ymin=lo, ymax=up)) +
	geom_hline(yintercept=0, color="white", size=1.5) +
	geom_ribbon(, alpha=.2) +
	geom_line(size= 1.75) +
	scale_x_continuous(limits=c(-2,7),  breaks=c(-2,-1,0,1,2,3,4,5,6,7), minor_breaks=NULL) +
	scale_y_continuous(limits=c(-.25,.75)) +
	xlab("Change in GDP") + ylab("Difference in probability of in- and out-partisans saying conditions got better") +
		ggtitle("(a) by economic conditions") +
		theme(panel.background=element_rect(fill="gray91"), 
		axis.text=element_text(size=10, color="black") ,
		axis.text.y=element_text(size=8, color="black") ,
		plot.title=element_text(hjust=.5, face="bold"))
g1 +	geom_rug(aes(x=econ), data=e.dat, inherit.aes=F, size=.5, alpha=.6)
ggsave("Fig3_a.pdf", width=5.5)







# FIGURE 3B

myX <- matrix(0, ncol=nrow(fixef(MC1)), nrow=1)
colnames(myX) <- row.names(fixef(MC1))
myX <- as.data.frame(myX)
myX$"(Intercept)" <- 1
myX$age <- 4.52
myX$white <- 1
myX$income <- 3
myX$educ <- 3
myX$female <- 1
myX$pres_election <- 1

myX1 <- myX
myX1$pidiOutParty <- 1
myX1$"pidiOutParty:pres_election" <- 1

fd.polar.MC1 <- matrix(NA, nrow=length(p.range), ncol=3)

for(i in 1:length(p.range)){
	OutX <- myX1
	OutX$polar_congress.m <- p.range[i]
	OutX$"pidiOutParty:polar_congress.m" <- p.range[i]
	OutX <- as.matrix(OutX)
	
	InX <- myX
	InX$polar_congress.m <- p.range[i]
	InX <- as.matrix(InX)

	ans <- fd.MCMCglmm.3cat(MC1, InX, OutX, Z=NULL, use="all", type="response")
	fd.polar.MC1[i,] <- as.matrix(ests.MCMCglmm.3cat(ans)[3,])
	
	print(i/length(p.range))

}

fd.polar.MC1 <- cbind(p.range.orig, fd.polar.MC1)
fd.polar.MC1 <- as.data.frame(fd.polar.MC1)
names(fd.polar.MC1) <- c("p.range.orig", "fd", "lo", "up")

g1 <- ggplot(data=fd.polar.MC1, aes(y=fd, x= p.range.orig, ymin=lo, ymax=up)) +
	geom_hline(yintercept=0, color="white", size=1.5) +
	geom_ribbon(, alpha=.2) +
	geom_line(size= 1.75) +
	scale_x_continuous(limits=c(.5,.75), breaks=c(.5,.55,.6,.65,.7,.75,.8,.85,.9), labels=c(".50",".55",".60",".65",".70",".75",".80",".85",".90"), minor_breaks=NULL) +
	scale_y_continuous(limits=c(-.25,.75)) +
	xlab("Elite polarization") + ylab("Difference in probability of in- and out-partisans saying conditions got better") +
		ggtitle("(b) by elite polarization") +
		theme(panel.background=element_rect(fill="gray91"), 
		axis.text=element_text(size=10, color="black") ,
		axis.text.y=element_text(size=8, color="black") ,
		plot.title=element_text(hjust=.5, face="bold"))
g1 +	geom_rug(aes(x=polar), data=p.dat, inherit.aes=F, size=.5, alpha=.6)
		
ggsave("Fig3_b.pdf", width=5.5)

round(fd.polar.MC1[1,],2)
round(fd.polar.MC1[nrow(fd.polar.MC1),],2)





# FIGURE 3C

InLo <- matrix(0, ncol=nrow(fixef(MC3)), nrow=1)
colnames(InLo) <- row.names(fixef(MC3))
InLo <- as.data.frame(InLo)
InLo$"(Intercept)" <- 1
InLo$age <- 4.52
InLo$white <- 1
InLo$income <- 3
InLo$educ <- 3
InLo$female <- 1
InLo$pres_election <- 1

OutLo <- InLo
OutLo$pidiOutParty <- 1
OutLo$"pidiOutParty:pres_election" <- 1


InHi <- InLo
InHi$awareness <- 1
InHi$"awareness:pres_election" <- 1


OutHi <- InHi
OutHi$pidiOutParty <- 1
OutHi$"pidiOutParty:awareness" <- 1
OutHi$"pidiOutParty:pres_election" <- 1
OutHi$"pidiOutParty:awareness:pres_election" <- 1

fd.polar.aware.MC3.lo <- matrix(NA, nrow=length(p.range), ncol=3)
fd.polar.aware.MC3.hi <- matrix(NA, nrow=length(p.range), ncol=3)


for(i in 1:length(p.range)){
	X.InLo <- InLo
	X.InLo$polar_congress.m <- p.range[i]
	X.InLo <- as.matrix(X.InLo)

	X.OutLo <- OutLo
	X.OutLo$polar_congress.m <- p.range[i]
	X.OutLo$"pidiOutParty:polar_congress.m" <- p.range[i]
	X.OutLo <- as.matrix(X.OutLo)

	X.InHi <- InHi
	X.InHi$polar_congress.m <- p.range[i]
	X.InHi$"awareness:polar_congress.m" <- p.range[i]	
	X.InHi <- as.matrix(X.InHi)

	X.OutHi <- OutHi
	X.OutHi$polar_congress.m <- p.range[i]
	X.OutHi$"awareness:polar_congress.m" <- p.range[i]	
	X.OutHi$"pidiOutParty:polar_congress.m" <- p.range[i]
	X.OutHi$"pidiOutParty:awareness:polar_congress.m" <- p.range[i]	
	X.OutHi <- as.matrix(X.OutHi)


	ans.lo <- fd.MCMCglmm.3cat(MC3, X.InLo, X.OutLo, Z=NULL, use="all", type="response")
	fd.polar.aware.MC3.lo[i,] <- as.matrix(ests.MCMCglmm.3cat(ans.lo)[3,])
	ans.hi <- fd.MCMCglmm.3cat(MC3, X.InHi, X.OutHi, Z=NULL, use="all", type="response")
	fd.polar.aware.MC3.hi[i,] <- as.matrix(ests.MCMCglmm.3cat(ans.hi)[3,])

	print(i/length(p.range))

}

fd.polar.aware.MC3.lo <- cbind(p.range, p.range.orig, fd.polar.aware.MC3.lo)
fd.polar.aware.MC3.hi <- cbind(p.range, p.range.orig, fd.polar.aware.MC3.hi)
fd.polar.aware.MC3 <- rbind(fd.polar.aware.MC3.lo, fd.polar.aware.MC3.hi)
fd.polar.aware.MC3 <- as.data.frame(fd.polar.aware.MC3)
names(fd.polar.aware.MC3) <- c("p.range", "p.range.orig", "fd", "lo", "up")

fd.polar.aware.MC3$awareness <- c(rep("Lowest     ", length(p.range)), rep("Highest", length(p.range)))
	
g1 <- ggplot(data=fd.polar.aware.MC3, aes(y=fd, x= p.range.orig, ymin=lo, ymax=up, linetype= awareness, group= awareness)) +
	geom_hline(yintercept=0, color="white", size=1.5) +
	geom_ribbon(, alpha=.2) +
	geom_line(size=1.75) +
	scale_x_continuous(limits=c(.5,.75), breaks=c(.5,.55,.6,.65,.7,.75,.8,.85,.9), labels=c(".50",".55",".60",".65",".70",".75",".80",".85",".90"), minor_breaks=NULL) +
	scale_y_continuous(limits=c(-.25,.75), breaks=c(-.25,0,.25,.5,.75)) +
	scale_linetype_manual(values=c(1,6), breaks=c("Lowest     ", "Highest"), name="Political awareness") +
	xlab("Elite polarization") + ylab("Difference in probability of in- and out-partisans saying conditions got better") +
		ggtitle("(c) by elite polarization and political awareness") +
		theme(panel.background=element_rect(fill="gray91"), 
		axis.text=element_text(size=10, color="black") ,
		axis.text.y=element_text(size=8, color="black") ,
		plot.title=element_text(hjust=.5, face="bold"))
		
g1 + 	geom_rug(aes(x=polar), data=p.dat, inherit.aes=F, size=.5, alpha=.6) +
		theme(legend.position=c(.17,.92), 
			legend.direction="vertical", 
			legend.key.width=unit(1.4, "cm"), 
			legend.text=element_text(size=10), 
			legend.title=element_text(size=10)) + 
		guides(linetype=guide_legend(title.position="top", title.hjust=.5, title="Political awareness", override.aes=list(size=1, fill=NA))) 
ggsave("Fig3_c.pdf", width=5.5)
dev.off()